Available online at www.sciencedirect.com 


ScienceDirect 


JOURNAL OF 


www.elsevier.com/ locate /jpowsour 


ELSEVIER 


Journal of Power Sources 171 (2007) 585—600 


Two-dimensional dynamic simulation of a direct internal 
reforming solid oxide fuel cell 
Jun Li *, Guang-Yi Cao, Xin-Jian Zhu, Heng-Yong Tu 


Institute of Fuel Cell, Department of Automation, Shanghai Jiao Tong University, Shanghai 200030, PR China 


Received 17 May 2007; received in revised form 10 July 2007; accepted 10 July 2007 
Available online 19 July 2007 


Abstract 


This study presents a two-dimensional mathematical model of a direct internal reforming solid oxide fuel cell (DIR-SOFC) stack which is based 
on the reforming reaction kinetics, electrochemical model and principles of mass and heat transfer. To stimulate the model and investigate the 
steady and dynamic performances of the DIR-SOFC stack, we employ a computational approach and several cases are used including standard 
conditions, and step changes in fuel flow rate, air flow rate and stack voltage. The temperature distribution, current density distribution, gas species 
molar fraction distributions and dynamic simulation for a cross-flow DIR-SOFC are presented and discussed. The results show that the dynamic 
responses are different at each point in the stack. The temperature gradients as well as the current density gradients are large in the stack, which 
should be considered when designing a stack. Further, a moderate increase in the fuel flow rate improves the performances of the stack. A decrease 
in the air flow rate can raise the stack temperature and increase fuel and oxygen utilizations. An increased output voltage reduces the current density 


and gas utilizations, resulting in a decrease in the temperature. 
© 2007 Elsevier B.V. All rights reserved. 
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1. Introduction 


The solid oxide fuel cell (SOFC) functions at high tempera- 
ture viaa completely solid-state ion-conducting electrolyte [1,2]. 
With the capability of internal reforming, the SOFC can be 
directly fueled with pure hydrogen, natural gas, coal gas and 
other hydrocarbons [1—3]. Further, the exhaust gas generated 
by the SOFC contains a lot of heat, steam and unused fuel, 
which can be conveniently recycled in the power plant when 
the SOFC is in use with gas turbines. Therefore, with such high 
energy generating efficiency and low pollution, the SOFC has 
been considered to be one of the most promising fuel cells for 
future power plants [1,2]. 

Several models have been proposed to simulate and analyze 
the operating conditions of SOFC stacks [4—15]. Costamagna 
and Honegger [8] presented and validated a simulation model 
for a SOFC stack integrated with an air pre-heater. Padulles et 
al. [9] demonstrated two other SOFC stack models to simu- 
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late the power system and the power conditioning units of the 
plant, respectively. However, dynamic responses of the temper- 
ature and current density in these simulated models are not clear. 
Recently, a novel one-dimensional dynamic model of a tubular 
SOFC with external and internal reforming was proposed by 
Jiang et al. [10]. Using this dynamic model, SOFC character- 
istics such as the cell voltage, temperature and power under 
different conditions can be predicted. However, the distribu- 
tion of parameters such as current density and gas concentration 
across the cell is not known. 

Compared with indirect internal reforming SOFC (IIR- 
SOFC) and hydrogen-fueled SOFC stacks, a direct internal 
reforming solid oxide fuel cell (DIR-SOFC) stack has unique 
reaction kinetics, in which the reforming reactions interact with 
concurrent electrochemical reactions [14—18]. It is well known 
that the dynamic performances and the parameter distribution of 
the stack are crucial for designing and controlling DIR-SOFC. In 
this regard, a one-dimensional DIR-SOFC stack model (anode- 
supported with intermediate temperature) has been used for 
simulation of DIR-SOFC [14,15]. Further, using this model, 
the one-dimensional steady-state parameter distribution and 
dynamic response to several current density step-changes can 
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Nomenclature 
A stoichiometric matrix 
Ak, pre-exponential factor for rate coefficient kn 


Aga pre-exponential factor for К ad 

B coefficient matrix 

Cp specific heat capacity (КЈ Ке! K^! for solid, 
kJ тої! K^! for gas) 


С p specific heat capacity vector 

E Nernst voltage (V) 

E, activation energy of reaction n (kJ mol`!) 

Fear Faraday constant (96,485 C шо!) 

Ё gas flow rate vector 

Ag" change of Gibbs free energy (kJ то! !) 

h gas enthalpy (kJ то! !) 

АН enthalpy change of adsorption (kJ тої!) 

i current density (A m7?) 

ig exchange current density (А m^?) 

iL limiting current density (А m~?) 

J gas component 

Kn rate coefficient for reforming reaction n 

K equilibrium constant 

K at adsorption constants for gas component j 

m number of moles of electrons transferred per mole 
of reactant 

n reaction 

P partial pressure (bar) 

q heat generation or transfers (kJ m^? s7!) 

r reaction rate (mol s~! m^?) 

T reaction rate vector 

R gas constant (8.314 J K^! mol`!) 

Sh volumetric rate of heat generation (kJ m`? s7!) 

T temperature (K) 

T temperature vector 

V output voltage (V) 

v gas velocity vector 

Greek symbols 

a*", а charge transfer coefficients 


Bi, fo coefficients for 7°™ 


у", у coefficients for i^ and ip^ 

material thickness (m) 

thickness vector 

emissivity 

thermal conductivity (W m`! K7!) 

ohmic polarization (V) 

activation polarization (V) 

concentration loss (V) 

convective heat transfer coefficient (W m~? K^!) 
density (kg m^? for solid, mol m^? for gas) 
density vector 

Stefan—Boltzmann constant (5.676 x 

1078 W m^? К“) 


З З I3 vw о ооо 


DID nmm 


3 


Superscripts 

an anode 

ca cathode 

in inlet 

Ip lower bipolar plate 

out outlet 

pa upper bipolar plate surface 
pe lower bipolar plate surface 
re reforming 

S SOLID (anode-electrolyte-cathode) 
sa anodic surface of SOLID 
sc cathodic surface of SOLID 


up upper bipolar plate 


Subscripts 

I reaction I 

II reaction II 

III reaction III 

conv convection 

e electrochemical reaction 
rad radiation 


be generated. However, because the two-dimensional situation 
is not considered, the one-dimensional model is only applicable 
to co- and counter-flow operation but not to cross-flow operation. 

Inthis study, we have constructed a two-dimensional dynamic 
model of the DIR-SOFC. This dynamic model is generated 
based on the methane steam-reforming kinetics, mechanism of 
chemical-to-electrical conversion, energy balance and mass con- 
servation. To describe the dynamic behavior in the cross-flow 
stack of this model, we have utilized multiple partial differen- 
tial equations (PDEs). We have also employed an alternative 
computational method to find the numerical solutions for the 
PDEs. Further, several cases are presented to study the steady 
and dynamic performances of the DIR-SOFC using this two- 
dimensional dynamic model. 

The present study is organized into five sections. In Sec- 
tion 2, the mathematical model is established following a brief 
description of the DIR-SOFC. In Section 3, a numerical method 
is used to solve the PDEs obtained in Section 2. In Section 4, the 
steady-state distributions (the temperature, current density and 
gas species molar fractions), as well as the dynamic responses 
of the stack for four cases, are presented and discussed. Sec- 
tion 5 concludes the paper and gives a perspective on the future 
research. 


2. Modeling of a DIR-SOFC 


The working principles of a DIR-SOFC stack are illustrated 
in Fig. 1. The electrolyte, which divides the stack into two elec- 
trodes, acts as an electronic barrier and avoids the direct chemical 
reaction of the fuel at the anode with the oxygen at the cathode. 
At the cathode, molecular oxygen combines with electrons and 
is reduced to negatively charged ions (O?-) with the aid of a 
catalyst [1,2]. 
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Bipolar plate 


Anode 
Electrolyte 


Cathode 


Bipolar plate 


Fig. 1. Diagram of a DIR-SOFC stack. 


Cathodic electrochemical reaction: 
1 _ 2— 
392 +2e — О (1) 


The steam—methane gas mixture flows in the anode of the DIR- 
SOFC and be reformed inside to provide the hydrogen and 
carbon monoxide for the electrochemical reaction to generate 
electricity [1,2]. The reactions at the anode of the DIR-SOFC 
are more complex than that using the pure hydrogen fuel due 
to coupling between reforming reactions and electrochemical 
reactions. 

At the anode, hydrogen and carbon monoxide react with oxy- 
gen ions to form water and carbon dioxide, respectively, and 
release electrons. 

Anodic electrochemical reaction: 


Н, + 027 > H50 + 2e7 (2) 
CO + O?- > CO; + 2e7 (3) 


Electrons released at the anode pass through the bipolar plate to 
the cathode of the next cell in the series, or follow the external 
circuit to drive the load, before finally returning to the cathode. 

Bipolar plates, with machined flow channels, clamp together 
the anode-electrolyte-cathode construction and provide the 
main structural support of the fuel cell. They also collect the 
current generated by the electrochemical reaction. 

Some assumptions are necessary in the derivation of the two- 
dimensional DIR-SOFC model. 


(1) All gases are assumed to be ideal gases. 

(2) The anode, electrolyte and cathode layers, which are 
clamped together, are assumed to be horizontally homo- 
geneous. The anode-electrolyte-cathode is denoted by the 
abbreviation SOLID. 

(3) The walls of gas channels are smooth. 

(4) Constant pressure is maintained within the gas channels. 


(5) Methane does not directly participate in the electrochemical 
reaction at the anode. 
(6) The stack works in the cross-flow mode. 


2.1. Reforming model 


The principal reactions in the methane steam-reforming are 
listed as follows [16—18]: 


Reaction I: CH4 + 2H50 <> CO» + 4H2 (4) 
Reaction II: CH4 + H20 <> CO + ЗН» (5) 
Reaction Ш: CO + H20 <> CO2 + H2 (6) 


In the reforming process, reactions I and II are the dominant 
reactions with products mainly comprising H2, CO and СО». 
Both reactions are endothermic. Reaction III is the water-gas 
shift reaction. The electrochemical reaction is exothermic and 
supplies heat and steam to the endothermic reactions. Below 
900 K, Reaction I dominates and the methane is mostly con- 
verted to carbon dioxide. Above 900 K, reaction II becomes 
predominant. The reversible reactions obey the laws of chem- 
ical equilibrium. The equilibrium constants for reactions I-III 
can be computed by three functions of temperature [11] which 
are given by 


Pco, P4 22, 430 
к= а ёк (- 4 26.078) | (7) 
Рсн, Pho T 
Pco P3 26, 830 
Kur ae (- + 30.114) ; (8) 
Рсн, Рн,о T 
Pco, Р, 4400 
Km = ‚СОН; Ho = ехр (+ = 4.036) А (9) 
РсоРн,о Т 


where Pcu,, Рн,, Рсо, Pco, and Рн,о are the partial pres- 
sures of СНД, H5, CO, CO» апа H20, respectively, when the 
chemical reactions reach equilibriums. T' is the temperature. As 
can be seen from Eqs. (7)-(9), the consumption of the hydrogen 
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caused by the electrochemical reaction can shift the reforming 
reaction in the forward direction. The methane steam-reforming 
is preferably carried out under the conditions of high tempera- 
ture and high steam-to-carbon ratio (S/C). Nickel in the SOFC 
stack is used as the internal reforming catalyst to accelerate the 
reforming reactions. 

The reforming reaction rates are computed by following 
equations proposed in [16]. 


k | Е рап pan 

P= | PP |, (10) 
Pip (DEN) Кү 
and 

re = EN. M pan pan _ Pity Poo (11) 
T i^^ (DEN)? „л, Ky | 
nis RS km an an __ Pip Poo, (12) 
Вавр \ S0 O Km /- 


The superscript “ап” denotes the anode. k,, with n = LILIII, 
are the rate coefficients for the reforming reactions and are 
calculated by 


E 
kn = Ap, ехр| ——— |], n= LILIII, (13) 
n р RT 


where Ар, is the pre-exponential factor of the rate coefficient 
kn, En is the activation energy of the reaction, and R is the gas 
constant (8.314 J K7! mol~!). DEN is defined as 


DEN = 1 + Ka Po + Ki, Pip + Kcu, Pen, 


Кз Ро 14 
7 H20 pan (14) 
H2 


In Eq. (14), кг, Kit. Ke. and Kilo. the adsorption con- 
stants for CO, H5, CH4and H20, respectively, are computed 
by 


а АН 
К? = Арка exp | — : , j=CO,H2, СНА, H20, 
. j 


RT 
(15) 


where A каа is the pre-exponential factor of the adsorption con- 


J 
stant Ка, апа АНА is the enthalpy change of adsorption. 
2.2. Chemical-to-electrical conversion of SOFC 


The electrochemical oxidation of CO is not considered in our 
study, because data in the literature show that the electrochemical 
oxidation of hydrogen is faster than that of CO [18]. Hence, the 
Nernst voltage Е is mainly dependent on the electrochemical 
reaction of H2: 


св RT (PRR? 
E = + n an , 
2 Fear Ро 


(16) 
2 Ера 


where Ag. is the change in molar Gibbs free energy of for- 
mation at standard pressure, and Ffar is the Faraday constant 
(96,485 C mol [1]. The superscript "ca" denotes the cathode. 
The pressure unit is bar. 


The cell output voltage, reduced by the ohmic polarization 


1°, the activation polarization 1! and the concentration loss 
n°", can be written as 
V=E- go E pe = ее. (17) 


The ohmic polarization can be expressed as 


gm = i Q, (18) 


Q = 8. B, exp (2). (19) 
where 6 is the material thickness, and i is the current density. 
The coefficients 6; and 62 can be obtained by curve fitting the 
simulated data proposed in Ref. [8]. The ohmic polarization 
depends on the physical configuration of the SOFC and on the 
operating temperature. 

The electrode activation polarizations can be described by 
the Butler—Volmer equation [8]: 


- an Far | act ca FFar ас 
i= ip len (« RT" ) exp (« RTs” ) Р (20) 
where o?" апа a“ are the charge transfer coefficients of the 
anode and cathode, respectively. The superscript "s"denotes the 
SOLID. The exchange current density io, which reflects the elec- 
trochemical reaction rate at the equilibrium potential, should be 
considered at both anode and cathode. 

Eqs. (21) and (22), proposed by Yamamura and Mogensen 
et al. [8,12,13,19,20], are used to calculate the exchange current 
densities at the anode and cathode, respectively: 


р рап рап pan 

п = у" ( ) ( чо) ехр (- =) : (21) 
Pret Pret RTS 
рса 0.25 Ee 

ig = y” 8 exp | ——& |, (22) 
Pref RTS 

where y?" and у“ are coefficients, Ей, and Е, are the activa- 


tion energies. P,ef is the reference pressure which is always set 
as 1 bar. 

As can be seen from Eqs. (21) and (22), the temperature and 
pressure influence the exchange current density. The electro- 
chemical reaction can be accelerated by increasing the operating 
temperature and the partial pressures. 

A formula related to the temperature and the current density 
is used to estimate the concentration loss according to the Nernst 
equation [2]: 


тое = LUE n (1 - ) (23) 
mF Far iL 
where ij, is the limiting current density. 

The consumption rates of Нә and О» due to the electrochem- 
ical reaction (rjj, ro.) are determined by the current density. 
The consumption rate of Н» is twice that of O2. A larger current 
density results in faster consumption rates of H2 and Ол in the 
electrochemical reaction. 


an __ === (24) 
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i 
i 25 
(027 двр _ 


2.3. Mass conservation 


The gas flow rates at the anode and cathode will change with 
the time and flow distance due to fuel and oxygen consumptions 
in the reforming and electrochemical reactions. 


ЕК" 4 A.7-B= po (26) 


F is the gas flow rate vector including the vectors of the anodic 
and cathodic flow rates. 


MEI 

P= |.|. (27) 
where 

> T 

Е = (Есн,, Fco, Fed,» Fipo Ён) ; (28) 
> T 

Е = (Fo,, Fry, Fino: Feon) - (29) 


7 is the reaction rate vector including the rates of reactions at the 
anode and cathode. 


yan 
? = n (30) 
where 
puce HE. (31) 
г = го. (32) 


А is the stoichiometric matrix. 


А? 0 
А = g awl (33) 
where 
=1 -1 0 0 
0 1 -1 0 
A™ = | 1 0 1 0 |, (34) 
—2 -1 -1 1 
4 3 1 -1 
—1 
0 
Аба = "s (35) 
0 


Bis the coefficient matrix in which elements are functions of the 
reaction area. 


2.4. Energy balance 
A considerable quantity of electrical energy is produced via 


electrochemical reactions. This is output through the polar plates 
and can supply electricity to external appliances: this is the 


motivation behind the development of fuel cells. As well as the 
electrical power, some energy is lost heating the stack and the 
rest flows out of the stack with the exhaust gases. 

The energy equation for the air and fuel flows can be 
expressed as 


Pah) + div (ov = 22) = Sh, (36) 
ot Cp 

where p is the density, Л the gas enthalpy, div denotes the diver- 
gence, V is the gas velocity vector, V denotes the gradient, ¢ is 
the thermal conductivity, C is the specific heat capacity and Sr 
is the volumetric heat generation source or sink. 

The energy equation for the SOLID and bipolar plate can 
be simplified as follows because of the incompressibility of the 
solid where there is no fluid motion: 


aT 
pCy- = div VT) = Sh. (37) 


The heat generation and transfers between the components 
of the stack are shown in Fig. 2. 

The local heat generated on the SOLID due to the reforming 
and electrochemical reactions can be expressed by 


Ф = -АН.- п Hp —AHgonp-Ved (38) 


where the subscript “e’’denotes the electrochemical reaction. 
The rate distributions of the endothermic and exothermic reac- 
tions are not uniform in the stack, so d? can be positive or negative 
at different positions in the stack. 

The heat change, caused by the water-gas shift reaction occur- 
ring in the anodic fuel flow, is 


ds = ~A Hm > тц. (39) 


The convective heat fluxes from the SOLID surface to the 
anodic and cathodic flows are calculated by 


ü = ё` (Т°*— Т“), (40) 
аа, = ЕТ Hn Т), (41) 


where £ is the convective heat transfer coefficient. Superscripts 
"sa"and "sc"denote the anodic and cathodic surfaces of the 
SOLID, respectively. 

The radiant heat fluxes between the SOLID and bipolar plate 
surfaces are calculated by 


9 4 
sup — op(T*” — ТҰ") 


Ягаа — ES a 1/euP " 1’ (42) 


slp _ ов(Т% = ТЇР?) 


Arad = туе ТЕР T _ 


where og is the Stefan—Boltzmann constant, and e is the emissiv- 
ity. Superscripts “up” and "Ip" denote the upper and lower bipolar 
plates, respectively. 

No reactions occur on the bipolar plates. Hence, we consider 
convective heat transfer, thermal radiation and heat conduction 
in the plates and neglect heat generation. 
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The convective heat fluxes from the upper and lower bipolar 
plate surfaces are given by 


üg = ёР°(Т"Р — Т“), (44) 
айка = gee(TIP — Та), (45) 


where superscripts “pa’’and “рс” denote the surfaces of the upper 
and lower bipolar plates, respectively. 

As Figs. 1 and 2 show, some heat (qfeact,n> qreaccn) 18 carried 
in the SOLID by the consumed reactant gases, and some heat 
(dřeact,n) 15 carried out of the SOLID by the product gases due 
to reactions at the interfaces of the SOLID and gas flow. The 
expressions of these heat fluxes corresponding to the various 
reactions are detailedly given in Appendix A. 


3. Computational method 


Problems in fluid dynamics usually require solutions of cer- 
tain PDEs with given boundary conditions. These equations can 
only be treated analytically in the case of simple boundaries. 
However, in the DIR-SOFC stack, the fluid governing equations 
with complex conditions are coupled and nonlinear. Therefore, 
it is difficult to acquire the exact mathematical solutions for 
these equations. Alternatively, a computational method is used 
to find an approximate numerical solution. First, the continuous 
physical domain is discretized into a set of cells to generate a 
computational grid. Second, the PDEs are transformed into finite 
difference equations. Lastly, the solutions of the difference equa- 
tions at grid nodes are used to approximate the variables of the 
PDEs [21]. 

Operating in the cross-flow mode, the fuel flows in the x 
direction, while the air flows in the y direction which is per- 
pendicular to the x direction. A regular rectangular grid, with 
grid spacing Au in the direction of x and Aw in the direction 
of y, is generated. A regular grid simplifies the discretization of 


Bipolar plate 


Anodic flow an 
d Shift 


s,an 
1 сопу 


the PDEs, saves computing time, reduces storage requirements, 
and improves computational accuracy. Combining the formulas 
derived in Section 2, the model equations are written in matrix 
form to facilitate programming and processing in the computer 
(see Appendix A for details). 


4. Simulation and discussion 


The two-dimensional dynamic model is developed in Mat- 
lab. Four cases are used to obtain the steady-state parameter 
distributions and the dynamic responses of the cross-flow DIR- 
SOFC under different conditions. Points with dimensionless (x, 
y) co-ordinates (0.025, 0.05), (1, 0.05), (0.5, 0.5), (0.025, 1), 
and (1, 1) are chosen as examples to show the dynamic per- 
formances at different positions on the SOLID. In addition, the 
mean temperature, the mean current density and the fuel and 
oxygen utilizations are important characteristics which must be 
determined. 


Case 1: Using the equations derived in Section 2 and Appendix 
A, the iteration is started with the initial values listed in 
Table 1. During the simulation, conditions in the stack 
evolve after each iteration until they reach a steady state. 

The two-dimensional steady-state temperature dis- 
tribution on the SOLID is illustrated in Fig. 3, and 
the corresponding contours are shown in Fig. 4. The 
two-dimensional current density distribution and the 
corresponding contours are displayed in Figs. 5 and 6, 
respectively. The molar fraction distributions of CH4, 
CO, CO», H20, Hzand О» are shown in Figs. 7-12 , 
respectively. 

The methane steam-reforming rate is rapid near the 
fuel inlet because of the large methane concentration 
(Fig. 7) and the function of the catalyst. The electro- 
chemical reaction rate of the hydrogen is slow at the 


up,an 


conv 


s 
| q react 


s,up 
rad 


q 


SOLID q* 


an 
q react," 


ca 
| q геасіл 


Cathodic flow 


q Ip,ca 


conv 


Bipolar plate 


| 


Fig. 2. Diagram of heat generation and transfers іп the stack. 


Table 1 
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Initial values of variables and parameters 


Item 


Value 


Fuel flow rate 
Air flow rate 


1.68 x 107? mol s~! 
1.43 x 1072 mol s^! 


Inlet fuel temperature 1023K 
Inlet air temperature 1023K 
Upper bipolar plate temperature 1024K 
Lower bipolar plate temperature 1024K 
SOLID temperature 1024K 
Output voltage 0.7V 
Pressure 3bar 
Inlet molar fractions at the anode 
CH4 28% 
CO 196 
CO» 196 
H20 60% 
H2 10% 
Inlet molar fractions at the cathode 
О 23% 
№ 75% 
CO» 196 
H20 1% 
inlet due to the low hydrogen concentration (Fig. 11). 
As a result, the local current density near the fuel inlet 
is very small (no more than 500 A т: Fig. 6). The 
heat generated in the electrochemical reaction is insuf- 
ficient to compensate for the thermal losses caused by 
the endothermic reforming reactions. This results in the 
temperature drop of the SOLID near the fuel inlet as 
shown in Figs. 3 and 4. The minimum local tempera- 
ture is 998 K which is lower than the inlet temperature 
(1023 K) of the fuel. 

The continuous reforming increases the hydrogen 
concentration along the fuel flow direction (Fig. 11), 
which results in the acceleration of the electrochemi- 
cal reaction as well as the rapid increase of the local 
current density. The reforming rate slows down due 
to the decrease of methane concentration as shown in 
Fig. 7. Then, the heat generated by the electrochemical 
reaction exceeds the reforming thermal losses. Thus, 
the temperature gradually rises and finally surpasses 

1300 + =, 
= 
e 


Air flow direction 


Fuel flow direction 


Fig. 3. Two-dimensional temperature distribution in the DIR-SOFC stack. 


Air flow direction 


i (A m^?) 


Fig. 5. 


Air flow direction 


Fuel flow direction 


Fig. 4. Temperature (K) contours corresponding to Fig. 3. 


0.4 


0.6 
0.8 


Fuel flow direction Air flow direction 


Two-dimensional current density distribution in the DIR-SOFC stack. 


the inlet flow temperature along the fuel flow direction 
(Figs. 3 and 4). The raised temperature also promotes 
the electrochemical reaction to increase the current den- 
sity and to generate more heat. 


Á LL ier: n — 4 - = 
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 
————- 


Fuel flow direction 


Fig. 6. Current density (A m^?) contours corresponding to Fig. 5. 
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Fig. 10. H2O molar fraction distribution. 
Most of the methane has been reformed, and the 
hydrogen concentration begins to decrease, once the air outlet (Figs. 5 and 6). The local temperature also 
fuel reaches the second half of the fuel channel declines from the maximum (1257 K) along the same 
(Figs. 7 and 11, respectively). Therefore, the local cur- direction, as shown in Figs. 3 and 4. 
rent density gradually decreases from the maximum The temperature rises gradually along the direction 
(4668 Am~~)along the fuel flow direction near the of the air flow (Figs. 3 and 4), due to heat loss from 
0.2 
S 015 S 
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Fig. 8. CO molar fraction distribution. Fig. 11. Ho molar fraction distribution. 
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Fig. 9. СОз тоаг fraction distribution. Fig. 12. O» molar fraction distribution. 
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Fig. 13. Dynamic responses of temperature to the initial settings. 


the SOLID into the gas flow. This is one of the reasons 
why the maximum temperature is not at the point of the 
highest oxygen concentration, but near the air outlet, as 
shown in Figs. 3 and 12. 

The difference between the maximum and minimum 
temperatures (1257 and 998 K) is 259 K, which indi- 
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cates a strong temperature gradient in the SOLID. This 
gradient should be taken as an important factor in the 
design and control of the stack. 

Initial values in the computational algorithm are 
equivalent to the step changes used to stimulate the 
model. The iteration process is the response of the 
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Fig. 14. Dynamic responses of current density to the initial settings. 
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system to the step change. The dynamic temperature 
responses at the five representative points to the initial 
settings are plotted in Fig. 13. The dynamic current den- 
sity responses are shown in Fig. 14. As can be seen from 
Fig. 13, it takes about 20,000 s for the mean tempera- 
ture to reach the stable value (1110 К). However, the 
time duration from the beginning of the iterations to the 
steady state at each point is different. Conditions at the 
point (0.025, 0.05) near the fuel and air inlet approach 
the stable value rapidly, while other points take more 
time to stabilize. The plots of the current densities, in 
Fig. 14, are similar to those of the temperature except 
the point (1,1) where the current density rises at first, 
before falling back to a lower level. The explanation 
for the dynamic response at the point (1,1) is given 
below. 

At the beginning of the iterations, the mean tempera- 
ture of the SOLID is relatively low and the consumption 
rate of the hydrogen in the fuel channels also remains 
low. Consequently, a large quantity of unused hydro- 
gen reaches the fuel outlet, increasing the local current 
density near the outlet. As the mean temperature gradu- 
ally rises with time, the overall hydrogen consumption 
rate in the stack will increase. Therefore, the hydro- 
gen concentration at the point (1,1) becomes smaller, 
and the corresponding local current density gradually 
decreases following the initial increase, as shown in 
Fig. 14. 

The fuel flow rate increases by 15% relative to the orig- 
inal rate; the air flow rate is unchanged. The curve of 
the fuel flow rate is shown in Fig. 15. 
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Fig. 15. Curve of the fuel flow rate. 


The dynamic responses of the temperature and cur- 
rent density for Case 2 are plotted in Figs. 16 and 17, 
respectively. The mean temperature and mean cur- 
rent density increase from 1110K and 1701 A m^?to 
1115K and 1886 Am"? with the increase of the fuel 
flow rate. However, there are also changes to the tem- 
perature and current density curves near the fuel inlet 
and outlet. Near the fuel inlet (points (0.025, 0.05) 
and (0.025, 1)), more heat is removed by the accel- 
erated fuel flow. As a result the local temperature 
gradually decreases, leading to a deceleration of the 
electrochemical reaction, and the corresponding local 
current density decrease shown in Fig. 17. On the 
other hand, the hydrogen concentration as well as the 
local current density near the fuel outlet (points (1, 
0.05) and (1,1)) increases and the corresponding local 
temperature increases (Fig. 16) due to the accelerated 
electrochemical reaction. If the fuel flow rate is suffi- 
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Fig. 17. Dynamic responses of current density for Case 2. 
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Fig. 18. Dynamic response of fuel utilization for Case 2. 


ciently increased, the cooling effect will become the 
main reason for the reduction of the overall tempera- 
ture. 

Utilization curves for fuel and oxygen are shown 
in Figs. 18 and 19, respectively. Fuel utilization sud- 
denly decreases following the increase in fuel flow 
rate, from 67.95% to 64.83%, while the oxygen uti- 
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Fig. 19. Dynamic response of oxygen utilization for Case 2. 
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lization suddenly increases at this time from 20.8696 
to 22.8796. Both the fuel and oxygen utilizations then 
gradually increase because of the intensified mean cur- 
rent density, before stabilizing at 65.52% and 23.12%, 
respectively. The reduction of the steady-state utiliza- 
tion of the fuel does not imply a decrease in methane 
consumption. Conversely, more methane is reformed to 
hydrogen for oxidation in the electrochemical reaction, 
thus the oxygen utilization for Case 2 increases. 

The air flow rate decreases by 10% relative to the orig- 
inal rate, while the fuel flow rate is unchanged. Fig. 20 
shows the curve of the air flow rate. 

The dynamic responses of the temperature and cur- 
rent density for Case 3 are shown in Figs. 21 and 22, 
respectively. Near the fuel and air inlet, at the point 
(0.025, 0.05), the local temperature and current density 
decrease immediately after the step change in the air 
flow rate. The overall temperature in the stack gradu- 
ally increases due to the weakened cooling effect of the 
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Fig. 20. Curve of the air flow rate. 
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Fig. 21. Dynamic responses of temperature for Case 3. 


air flow. The mean temperature and mean current den- 
sity increase from 1110 K and 1701 A m^?to 1117.7 К 
and 1720.4 A тт, respectively. Fuel and oxygen uti- 
lizations increase, from 67.95% and 20.86% to 68.72% 
and 23.45% as shown in Figs. 23 and 24, respectively. 
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Because of the increased consumptions in the stack, 
the hydrogen and oxygen concentrations reduce at the 
point (1,1) near the fuel and air outlet, reducing the cor- 
responding local current density, as shown in Fig. 22. 
However, the corresponding local temperature at the 
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Fig. 22. Dynamic responses of current density for Case 3. 
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Fig. 24. Dynamic response of oxygen utilization for Case 3. 

point (1,1) (Fig. 21) initially increases due to the heat 

transfer, and then decreases due to the reduced rate of 

the local electrochemical reaction. 

Case 4: The output voltage changes from 0.7 to 0.8 V; the flow 


rates are unchanged. The output voltage is plotted in 
Fig. 25. 
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Fig. 25. Curve of the output voltage. 


Figs. 26 and 27 show the dynamic responses of the 
temperature and current density for Case 4, respec- 
tively. The temperature and current density generally 
decrease following the increase of the output volt- 
age. The mean temperature reduces from 1110 to 
993 K, and the mean current density from 1701 to 
623.5 А m 2. The heat generated is greatly reduced due 
to the decrease in mean current density. A substantial 
fraction of hydrogen and oxygen is not consumed in 
the electrochemical reaction, which increases the con- 
centrations of the hydrogen and oxygen near the outlet 
following the step change in output voltage. Conse- 
quently, the local current density at the point (1,1) 
increases for about 8000 s, and the resulting heat gen- 
erated raises the corresponding local temperature for 
around 3000 s. However, the overall drop of the stack 
temperature lowers the outlet temperature, decreasing 
the corresponding local current density. Fuel utilization 
decreases from 67.95% to 24.91% as shown in Fig. 28, 
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Fig. 26. Dynamic responses of temperature for Case 4. 
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Fig. 28. Dynamic response of fuel utilization for Case 4. 
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Fig. 29. Dynamic response of oxygen utilization for Case 4. 


and the oxygen utilization decreases from 20.86% to 
7.69% as shown in Fig. 29. 


5. Conclusions 


A two-dimensional dynamic DIR-SOFC model has been 
presented based on the reforming reaction kinetics, the elec- 


trochemical model and the transfer principles of mass and heat. 
Four cases, comprising standard conditions, a fuel flow rate step 
change, an air flow rate step change and a voltage step change, 
are used to stimulate the numerical model to obtain the steady 
states and dynamic responses of the DIR-SOFC. Distributions 
of temperature, current density and the molar fractions of chem- 
ical species; dynamic responses at representative points; and the 
variations in fuel and oxygen utilizations, are presented to reveal 
the performances of the cross-flow stack. The results show that 
the temperature gradients on the SOLID are sufficiently large 
that they should be considered in design and controlling of the 
stack. The fuel flow rate, air flow rate and voltage can influence 
the distributions of the stack parameters. The dynamic responses 
at different points on the SOLID are not the same. In Case 2, the 
increase in the fuel flow rate increased the mean current density 
and mean temperature; the fuel utilization decreased, while the 
oxygen utilization increased. However, when the fuel flow rate 
is sufficiently increased, the cooling effect will yield a decreased 
temperature. In Case 3, the mean temperature as well as the mean 
current density increased with the decrease in the air flow rate 
due to a weakened cooling effect. Therefore, the fuel and oxygen 
utilizations also increased. Results from these simulations show 
that the fuel and air flow rates should be balanced between the 
performance improvement and the temperature optimization in 
practical applications. In Case 4, the increase in the output volt- 
age decreased the mean current density and the generated heat. 
Hence, the mean temperature and utilizations both decreased. 

There are still many problems that have not been resolved 
completely in this area, such as the optimization and control of 
the DIR-SOFC. Future work will concentrate on the optimiza- 
tion of modeling and control, to improve the performances of 
the DIR-SOFC stack. 
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Appendix A. Matrix form representation of the model 
equations 


The model equations are presented in matrix form as follows: 


diag) - іар (л, у) · авс“, ,) - diag(T s.) 


= diag($1) + diag(»), (A.1) 
8 = (8°Р, 9n, 88, 8%, al)". (A.2) 
p = (о%?, p™, p, р, pP)", (A.3) 
(e qon qe (A.4) 
T = (їр, T*^, TS, 7%, TP)’, (A.5) 


where 8 is the thickness vector, p is the density vector, C p is the 
specific heat capacity vector, T is the temperature vector, diag(-) 
denotes the diagonal matrix and the subscripts “x, y"denote the 
co-ordinates (x, y). The matrices introduced on the right hand 
side of Eq. (A.1) are: 
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Foy 0 bs 
Тап 
аир (А.9) 
= 
o 17 1 
^u? —2 
0 TM; y 0 1 
5 ecs M , А.10 
91 Ü 0 | 0 TM; „1 1 ) 
Ашт? -2 
0 1 


Fea T rca 
& — Aulay! | 07 Созу 0 
9p = poa 0 Сса 


х,у р,х,у 


та _ 
shal (A.11) 
=, 
o 1” 1 
^u? —2 
Ip 
0 ТМ, 0 
g? = 8р = ЫТ ‚ (А12) 
0 0 TM, 1 
^w? —2 
0 1 
0 Ty, y-1 0 
ТМ, у == Ty-1,y Туу Ty , (A.13) 
0 Ty, y+1 0 
Фу. = (облу + Tad x. y (A.14) 
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05" = (eon + a Treacte)y,y* (A.17) 
I а ‚1 

Фр == (42,2 + rue y (А.18) 


In Eqs. (A.15)-(A.17), qreactn> бтеасьт 24 Фес Сап Бе 
written as a general equation to facilitate matrix operations: 


> > 5 2T 
= Dil. [diag(T) & C- C, )]- DER 


1 
dreact,n n^ 
n = LlLe, 


(A.19) 


1 = an,s,ca, 


where D^. and DI are the left and right factor matrixes, 
respectively, and & denotes Kronecker product. Each of these 
quantities has the following components: 


Cp = (С›,сн„. C pco; Cp.CO»; C p, H50» 


C pH»; C p.05; Comey s (A.20) 
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(1) The electrochemical reaction 
The components of D$. and D$ corresponding to dřeact,e 
are 


ài = ash = att = ast = (0,0,0,0, 0), 


(A.24) 


600 


(2) 


(3) 
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dst = (0,0, 0, 1,0), (А25) 
aR — gk — ДК = ae — 0,0,0,0,0,0,0), (A26) 
= (0, 0, 0, 1, 0, 0, 0). (A.27) 
The components of D2" and D2™® corresponding to deste 
are 
dm. = dm = qn = ganL = (0, 0, 0,0, 0), (A.28) 
al — (0, 0, 0, 1, 0), (A.29) 
дак _ дак = да — gen = (0, 0, 0, 0, 0, 0, 0), 
(A.30) 
Дау = (0, 0, 0, 0, 1, 0, 0). (A.31) 


ca,L ca,R Ё са 
The components of D£^* and D£^* corresponding to гелсе 
are 


derum = aoe = TR = (0, 0, 0, 0, 0), (A.32) 
dc^" = (0,0, 0,0, 1), (A.33) 
der = dd = аск = (0, 0, 0, 0, 0, 0, 0), 
(A.34) 
аса = (0, 0,0, 0, 0, 1,0). (A.35) 
Reaction I 
The components of Dj" and О corresponding to dieactI 
are 
а d; = (0,0,0,0,0), (A.36) 
dis = (1, 0, 0, 0, 0), (A.37) 
at = diy = Фр = dF = 0,0,0,0,0,0,0), (A38) 
= (0, 0, 1, 0, 4, 0, 0). (A.39) 
The components of ppt and Dj aR corresponding to esl 
are 
di = di3" = diy’ = dee" = (0,0,0,0,0), (A40) 
йу = (1, 0, 0, 0, 0), (A.41) 
d cdd ed = (0, 0, 0, 0, 0, 0, 0), 
(A.42) 
dio = (L, 0, 0, 2, 0, 0, 0). (A.43) 
Reaction II 
The components of Du and рК corresponding to 


s 
dreact,II are, 


ain = dy; = di4 = di 115 5 = (0, 0, 0, 0, 0), (A.44) 
a! = (0, 1,0,0,0), (A.45) 
aR — ask — aR — GR _ (0,0, 0,0, 0,0, 0), (A.46) 
dik = (0, 1, 0, 0, 3, 0, 0). (A.47) 


The components of D? and 2 corresponding to q% a 
are 


~an,L ~an,L qan L 


dr IL3 ^ "ILA ^ dy 5 = (0, 0, 0, 0, 0), (A.48) 
dy» = (0; 1,0,0,0); (A.49) 
йт = daa = dag = da3 = (0,0,0,0,0,0,0), (A50) 
dit’ = (1,0, 0, 1, 0, 0, 0). (A.51) 
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